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SUMMARY 


A numerical procedure has been developed to compute the inviscid super/ 
hypersonic flow field about complex vehicle geometries accurately and effi- 
ciently. A second-order accurate finite difference scheme is used to integrate 
the three-dimensional Euler equations in regions of continuous flow, while all 
shock waves are computed as discontinuities via the Rankine-Hugoniot jump 
conditions. Conformal mappings are used to develop a computational grid. The 
effects of blunt nose entropy layers are computed in detail. Real gas effects 
for equilibriiim air are included using curve fits of Mollier charts. 

Typical calculated results for shuttle orbiter, hypersonic transport and 
supersonic aircraft configurations are included to demonstrate the usefulness 
of this tool. 

A computer code utilizing this computational procedure is described in 
Volume II of this report. 
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INTRODUCTION 


The prediction of the steady three-dimensional inviscid flow fields about 
a vehicle is of great interest to the designer. Most data necessary to develop 
a high-speed vehicle is presently obtained from wind tunnel tests which are ex- 
pensive, slow and sometimes inadequate. The goal of this work was to create a 
computer code to be used in the development of supersonic vehicle configurations 
This code should therefore meet three basic requirements. The first is appli- 
cability. In order to obtain the required accuracy for the problem of computing 
the flow over a wide variety of geometries, for a wide range of Mach numbers and 
angles of attack, the Euler- equations must be solved (Fig. l). Small perturba- 
tion techniques yield accurate results only for the flow over slender bodies 
flying at low supersonic Mach numbers and small angles of attack, while 
Newtonian theory yields useful results only for large Mach numbers. Neither of 
these theories can yield all the details of the flow even in their range of 
applicability . 

The second requirement is efficiency. The calculation of the flow field 
over a complete vehicle should take no longer than two hours on the C.D.C. 6600 
computer. This requirement can only be met by reducing the number of mesh 
points needed to obtain an accurate solution and keeping the program logic as 
simple as possible. Computational techniques that "capture" the shocks in the 
flow field require too many mesh points to obtain acceptable results (ref. 1). 
The three-dimensional method of characteristics is rejected because of its 
extreme complexity in program logic. 

The last requirement is that the code should be a user-oriented tool. This 
is in contrast to codes that are tailor-made for a particular configuration 
(ref. 2) and codes that must be constantly monitored '‘to nurse the solution 
through critical regions" (ref. 3). The designer should only have to specify 
the vehicle geometry and flight conditions to obtain reliable results in a 
directly usable form (e.g., aerodynamic -coefficients , boundary layer inputs, 
etc.). The vehicle geometry should be input via techniques that, on the one 
hand, are of the same advanced level as those used in the best incompressible, 
supersonic and hypersonic three dimensional tools, and, on the other hand. 



possess longitudinal and cross sectional continuity needed when solving partial 
differential equations . 

Although the genei-al background of the numerics involved in solving this 
problem were available at the beginning of this study, no computer tool for 
actiially carrying out accurate flow field calculations past realistic con-' 
figurations existed. 

The only linnitation inherent in the present fomulation of this problem 
is that the Mach number in the marching direction (an axis running from the 
nose of the vehicle to its tail, figure 2) must be supersonic at every point in 
the flow field. This limitation implies that, first, the free stream Mach 
nixmber must be "sufficiently" supersonic and, second, that the geometry of the 
vehicle is such that there are no imbedded regions of subsonic flow. The region 
around the nose of - blunt nose vehicles can be computed using other existing com- 
putational techniques (e.g., ref. U) and once the flow becomes supersonic the 
present technique can he applied. In general this limitation means that' com- 
pressions in the flow field that cause subsonic Mach numbers cannot be handled. 
The general numerical scheme used to solve this problem has been developed by 
Moretti (refs. '5-8). ' It follows a number of basic guidelines: 

s A second order accurate finite difference marching technique (satisfying 
the C.F.L. stability condition) is used to numerically integrate the 
governing partial differential equations 

® xAll shock waves. in the flow field are followed and the Rankine-Hugoniot 
conditions are satisfied across them 

9 The intersection of two shocks of the same family is computed explicitly 

« Conformal mappings are used to develop a computational grid 

« The body boundary condition is satisfied by recasting the equations 
according to the concept of characteristics 

IS The edge of the entropy layer on blunt nose vehicles is followed from 
its origin and special devices are used to form derivatives across it 
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• Real gas effects are included (equilibrium air) when appropriate, 
by using curve fits of Mollier charts 

9 Sharp leading edge wings are computed using a local two-dimensional 
solution 

A computer code has been developed which uses these basic ideas. It is 
described in detail in Volume II of this report. This code has been used 
extensively to compute external flow fields and has been found to yield 
accurate results for a wide variety of complex vehicles (Fig. 2) flying at a 
wide range of supersonic Mach mmibers (M ^ 2 -* 26 ) and angles of attack 

09 

(o’ ~ + 30°)- Computed results are presented in this report to demonstrate 
and validate this Oomputational procedure. 

35 
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Figure 1. - Regions of applicability of inviscid flow theories 
(for the surface pressure on a sharp cone). 
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PROBLEM DEFINITION 


The problem of computing the external flow about a vehicle flying at 
supersonic speeds is completely defined by the vehicle geometry and the free 
stream Mach number, flow direction and ratio of specific heats (y) for an ideal 
gas (NOTE: throughout this report an ideal gas will be assumed; modifications 

of the computational technique for a real gas will be discussed in a later 
section) . 

In order to apply the present computational procedure, the free stream 
Mach number must be supersonic. It will be assumed first that the geometry has 
a plane of symmetry and second that the free stream velocity vector lies in this 
plane of symmetry. These two assumptions are not necessary but the present 
techniques have not" yet been applied to asymmetric bodies or bodies at yaw. 

In order to apply the boundary condition at the vehicle surface (i.e., 
vanishing of the nonnal velocity) the surface and all its first derivatives 
must be defined. Therefore an analytic definition of the geometry is needed, 
with continuous first derivatives. The second derivatives of the body geometry 
appear explicitly in the present formulation of the problem but continuity is 
not necessary. The flow field over geometries with discontinuous slopes, such 
as cone-cylinder combinations with sharp expansion corners or geometries with 
discontinuous wing-fuselage roots and/or canopy- fuselage connections (Fig. 3) 
have been computed with no special treatment. In these cases computational 
points on either side of the discontinuity have well defined slopes (the slope 
at the corner itself is defined by the limit from one side or the other) , 
therefore as far as the computation is concerned there is a smooth transition 
from one point to the next. Although major difficulties have not been en- 
countered in these cases, best results are obtained for vehicles on which these 
discontinuities are eliminated by fairings. Special treatment of other slope 
discontinuities such as sharp leading edge wings will be discussed later. 

Since a marching computational technique is employed (i.e., given data on 
a plane z = constant, (Fig. h), data on a plane z + A z is computed) a geometry 
definition that specifies cross sections in the x,y plane (Fig. 5 ) at each 

i is required. 


value of 



In order to define the body geometry as a single valued function of two 
variables, a polar coordinate system is used with the pole at a point in each 
cross section that will specify the body radius as a single valued function of 
the polar angle. Of course, the pole of this coordinate system is a function of 
the axial coordinate. A geometry package that fits all these requirements has 
been developed by A. Vachris and L. Yaeger (ref. 9). 

The main effort in this study is to develop a computational technique that 
marches downstream from an initial data plane (z = constant. Fig. 4) in which 
all the dependent variables are known and the component of the velocity in the 
z direction (w) is supersonic. 

For the flow near the nose of the vehicle there are efficient methods 
available so that the initial data can be easily computed. When the nose of 
the vehicle is blujit (Fig. 4a) the three-dimensional transonic blunt-body solu- 
tion developed by Moretti (ref. 4) can be used. If the nose is sharp (Fig. 4b) 
a region very close to the tip is conical. Input data from the flow field 
s(Dlution over sharp conical bodies with attached shocks are readily available 
(refs. 10 and ll). However, a method has been described by Moretti and Pandolfi 
(ref. 8) for computing sharp conical solutions. With this technique, a blunt 
cone solution will asymptote to the sharp cone solution, as the computation 
proceeds downstream. So that for sharp nose bodies, the calculation is started 
i^ith a blunt nose solution and continued downstream until the sharp cone solu- 
tion is reached and then the calculation is restarted with this sharp-nose 
solution. 




Figure 3- - Geometries with sharp corners. 
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Figure 4. - Coordinate system definition. 
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COMPUTATIONAL FRAME 


The definition of a computational grid is one of the most crucial steps 
in the building of a numerical technique. The idea of using conformal mappings 
to develop the computational mesh in this problem was originally proposed by 
Moretti (ref. 5). 

Three coordinate systems or spaces will be referred to: the physical 

space (x,y,z)j a mapped space (rj 0 ,p), and a computational space (X,Y,Z); 
figure 6. The physical coordinate system is Cartesian and defines the three 
velocity components that are computed. The governing equations are written 
in the physical space and then all derivatives are transformed into the com- 
putational space where the mesh points are at uniform intervals A^, A-% A^. 

Each constant plane in the mapped space is obtained by conformally 
mapping the geometry cross sections in the physical space into near circles, 
in the mapped space. So that the region bounded by the body and the bow shock 
becanes a ring in the mapped space. The corresponding computational space is 
obtained by normalizing the radial distance between body and shock (X-direction) 
and the circumferential distance between the two symmetry planes ( Y-direction) 
(Fig. 7). The mapped space serves three purposes: first, it distributes the body 
mesh points (which are evenly spaced in the computational plane) so that the 
necessary resolution is obtained in regions of large curvature where truncation 
'error may otherwise become too large. This is demonstrated in figures 7 and 8 
where the mesh points are concentrated near the tips of the wing and the tail 
in the physical space. Second, the mapping makes the body and bow shock posi- 
tions single valued functions r = b(6,^) and r = c(9,p) in the mapped plane. 
Third, since shock waves imbedded in the flow field become mesh lines, the only 
chance of success is that their traces in a mapped space are (for the crossflow 
shocks) predoiriinantly along lines 6 = constant and (for the wing-type shocks) 
predominantly along lines r = constant. This enables one to define a number of 
regions in the corajutatioiia.! space when more than one shock exists in a cross 
section. The dotted lines in figiue 9 are extensions of shocks to complete the 
boundaries of these regions. The points on either side of these portions of the 
boundaries are allowed to pass information across the boundary. The require- 
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merits for such a surface are that they meet the shocks at the edge points 
(Fig. 9) and that they do not intersect each other if they are of the same 
type. The cross flow surfaces (i.e., extensions of cross flow shocks) are 
defined as 9 = constant. The constant taken is the value of 9 at the last 
point on the shock. Since the wing-type shocks intersect each other, points 
on their dividing surfaces are taken at the same percent of the distance be- 
tween its two neighboring boundaries at its edge shock points, so that no 
point on the dividing surface of a wing-type shock gets any closer to the two 
adjacent boundaries of the shock than its edge shock points. 


Now the transformation (r,9,p) -» (x,y,z) is considered. A series of 
algebraic, invertible, conformal mappings have been developed which transform 
a wide variety of cross sections (Pig. 2) into near circles. 


Let Q = re^®, W. 


1,2, 3,4, 5, 


be intermediate spaces, 

, 

and G = X + iy (Fig. 10 ) 

.. i0 

(a) 

C - F/c 

(b) 

W^= + iE 

(c) 

W3= + B/W^ 

(d) 

W^= + iA 

(e) 

+ D^/(4Wi^) 

(f) 

G = + iC 

(g) 

0,^), y = while 

z* = ^ completes the 


definition of the transformation. 


Each coefficient A,B,C,D,E and F is determined by placing the singularities of 
the mappings inside the body so as to obtain a near circle in the g-plane. 
Therefore these coefficients are functions of the axial coordinate z. The 
coefficients are determined by the geometry. These mappings can handle a wide 
variety of vehicles with the coefficients defined as follows*. 
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A = 0 

(a) 

B = i(D3/54) 

(b) 

o 

II 

ro 

(c) 

B - (y^-y^)^ 

(d) 

E = (nii(W ) + Im(W^ )) /2 

(e) 

F = (W - iE)^/4 
5 

(f) 

where , W Wo the positions of points 1, 4 and 5 (Fig. lO) in 

h 2l, =5 

the Wg plane. These points are computed by inverting equations (l) as 

follows, ^ ^(y^ - y)2 _ 


and 



(a) 

=1, = iy4 

(b) ^3) 

Gj = iy. 

(c) 

and for K = 1,4 and 5 


"5^ = "k - 

(a) 

- 4w W. + = 0 

(b) (4) 

W = W, - iA 

(c) 

W ^ - W W? - B = 0 

(d) 

The roots with the largest modulus is used in equations 

(4b and d) . 

* A = 0 for all applications thus far, A ^ 0 would give 

the mappings 


even greater generality. 
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The derivatives A , B , etc., are obtained in a straight forward manner 
z z 

by differentiating equations (2-4) and are functions of , y^ , y^ , 

-,r^\ z ■^z z z 

and y^ (Fig. 10). 

z 

To obtain r,0 in terms of x,y(atp= z) all of the mappings must be in- 
verted. This is a straightforward process that can be carried out using 
equations (la) through (ig) (similar to the procedure used in equations (4)). 

In order to transform the derivatives in the governing equations from the 

physical space to the computational space the first and second derivatives of 

all the transformations are needed. The derivatives r , r , 0 , and 9 can be 

X y X y 

calculated as follows . 




and 


2 ""4 "5 

Ca = 1/G, 


=G - -Q 

If C = u + iv, where u = rcos0 and v = rsin6, then 

Ux - R6 (Cq) 

\ = MQq) 

and from the Cauchy-Riemann conditions 


(5) 

( 6 ) 


u = -V 

y X 


V = u 

y X ■ 


thus 


X 



w^)/r 

(a) 


(uu + 


(b) 


y 

y 


(T) 

(uvx - 

vux)/r 

(c) 


(uv - 

y 

vu )/r 

y ' 

(d) 
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For the ^-derivatives of the mappings the following results are obtained: 

= -Fp/C (a) 


= -Fp/C 

= w, „ 

+ iE„ 


P 

II 

^^3W2 

II 

oo^ 

+ iA 

II 

'M 

‘V, 


+ iC„ 

5P 

7 


y = Im (G ) (h) 

f r 

Where A„ = A , = B etc., since z = ^ and z = 1. Now since r,0, p are 

independent coordinates, it follows that ^ therefore 

= '■x yp * --z ° 

and t ^ ®y ^'p * \ "p = ° 

solving for r and 0 the following results are obtained 
z z 


r, = -(Xy y^ X x^) (a) 

(9) 

'’z = -<V yp * ''p) 

The second derivatives of this transformation needed for the body point 
calculation are derived in Appendix A. 
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The singularities of these mappings are all inside the body so that the 
mappings are never evaluated at singular points. It is not necessary that 
the m.appings he conformal but it has been found that conformal mappings give 
the best results in terms of m.esh point distribution. It is also not necessary 
that r, 9, p be orthogonal coordinates, in fact in general they are not. 

These mappings (equations (l))use simple algebraic expressions and their 
coefficients are defined explicitly so that transformation from one space to 
another takes a minimum of time. Considerable work has been done to develop 
conformal mappings that can map arbitrary cross sections into circles or near 
circles (ref. 12). These generalized mappings offer a greater flexibility than 
the mappings used herein but would require a large increase in ccmputational 
time. 

With the mapped plane completely defined, the transformation between the 
computational space and the mapped space (X,Y,z) (r,0,p) is required. Con- 

sider a cross section (Z = p= constant) with multiple shockwaves, e.g., 

(Fig. 9) bow shock, wing shock and tail shock plus two cross flow shocks. The 
computational plane is divided into IC regions in the Y direction, and LC 
regions in the X direction, they are ordered as in figure 9* body is 

described by r = B(Y,Z) and a wing-type shock as r = C (Y,Z) for = 1 -* LC 

A) 

(;J = LC being the bow shock). Similarly, the cross flow shocks are described 
as 9 = (X,z) for i = 1 -• IC + 1, (i = 1 is the bottom symmetry plane and 

i = IC + 1 is the top symmetry plane). These surfaces are shocks for some 
range of X and Y and arbitrary (dividing) surfaces for other values of X 
and Y as described previously. 


Wow define LC + 1 surfaces such that: 


C^(Y,Z) = B(Y,Z) 

^ ^ . . .LC + 1) 

transformation to the computational plane then can be written as: 
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X = - c^) 

(a) 


Y = (0-Hp/(H.^j^ - HP 

(b) 

(10) 

z = ; 

(c) 



The coordinates X,Y,Z are not orthogonal. The boundai'y C in the mapped 

X 

plane becomes X = 0 in region £, and becomes X = 1. Similarly, the IC 

regions in the Y-direction in the computational plane are bounded by 
Y = 0 and Y = 1. 


Inverting this transformation yields the result: 


r = x(c^_, -C ) + C 
i+1 i a 

(a) 


S = + H. 

(b) 

(11) 


= Z 


(c) 


Again the derivatives of this transformation X ,X ,X„,Y ,Y and Y are 


needed. 

First 


^x = (^,.1-^,) 

r^ = (C^ - C ) X + C 

i+1 i i 

^ 7 . = ■■ "z ) ^ "Z 

£+1 i i 

H - ( V ' * \ 

“ ^^i+1 ■ 

%= («z.,, - «z.) 

1+1 1 1 


r’^^e’^'-p’ r’ 9 " P 

(a) 

(b) 

(c) 


(e) 

(f) 


( 12 ) 


= 0 


P,= 0 


(g) 

(h) 


= 1 


(i) 
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The Jacobian of the transformation is defined as : 


j _ 5 . (^; 9}^ 
- a(x,Y,z 


®x ®y ®z 
^x ^z 


and thus 


X = 1 
J 


r ^ 


X. = 


Y. = 


0 T 


X, = 1 

J 


= 1 
? J 


1 

^ r\ 
1 

(a) 

_a(x,Y,z)_ 

(b) 

B(r,X,2) 

_a(x,Y,z)_ 

(c) 

_d(x,y;z)_ 

(d) 

d(r»9.X) 
_B(X,Y,Z) _ 

(e) 

5(r,e,Y) 

d(x,Y,z) 

(f) 


(13) 


After some algebraic manipulations the above derivatives can be written 
in the following form: 


X = 1./ 
r ' 


6 + XD^e^ + 


Y = D^X 
r 1 r 




*9 = 


(a) 

(b) 

(c) 

(d) 


(l4) 
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4 ° " 

’x 62 + X 6^3 + D 3 - 

(e) 


[6 + X 6 , + 

(14) 



(f) 

where 

^ " ^ji +1 ■ 

(a) 


= Si +1 - 

(b) 


" ^Zj^+l ■ ^Zi 

(c) 


A =«i.l-«i 

(d) 


■ ^X. 
1+1 1 

(e) 


^ \+l ^i 

(f) (15) 


Di = - (YAjj - Hx.)/a' 

(g) 


1^2 = - (X 6 y + 

(h) 


D - - (IL + Ya )/a 
0 ^i 

(i) 


= - <“x. * ^Ax>/a 

(j) 

This transformation 

is a modification of the one 

used previously by 


Moretti to so3.ve nizmerical problems, (e.g. see reference 4) Singularities 

occur when C -- C or H. = H. , ^ and when J = 0. The former occurs when 

i f 1 1 1 + 1 

two shocks intersect j this matter will be discussed in the section on 
"Treatment of Shocks". The latter case occurs when the mesh lines X = constant 
and Y = constant become parallel in the physical plane. This can occur for 
certain locations of cross flow shocks. However this problem can be overcome 
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by either modifying the conformal mappings so that the cross section in the 
mapped plane is "more circular" or using a cross flow shock t 3 ^e surface 
(which acts like an extension to a cross flow shock) to control the shape of 
the mesh lines. 

All shocks are defined in the mapped plane as r = c(0,2) and 0 = h(r,p) 

so that C(Y,Z) = c [e(X ,Y,Z),2], H(X,Z)=h [r(X,Y ,Z),2] (where X and Y 

s s * s s 

are either 0 or l) and their derivatives and 

The body is defined in the physical plane and an iterative procedure is 

needed to describe the body as r b(0,p) from which B(Y,Z) ca,n be computed. 

From the derivatives b , b , c , c_, h and h the calculation of B , C , B , 

0 2 0 2 2 r ^ 

and C proceeds as follows. 

Zi 

Using the notation of equations (ll - 13) define 

C^ = b(e,p) 


must be calculated. 



then 


^X. = ^r. 


- H(Xg, Z)^] 

C(Y^,Z)^] 


X+1 


(16) 


Where again X (the value of X at the shock C ) and Y (the value of 
s' Z s ^ 

Y at the shock H^) are 0 or 1. 






• 2 - 

1 




(a) 

1.1 

* Wh 

(b) 

^r. '■z 


(c) 

1 





(d) 


(IT) 
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At the points X = and Y = at a cross section Z these equations 
result in a set of simultaneous linear equations for IL,(X ,Z) and C„(Y ,Z). 

S ^6 S 






Ty (h 


’^i+l 


h ) - h ■ 
r. r. 

1 i-j 


(18) 


where as H^Cx^jZ) can be computed from equations (l6). For all other points 
equations (l6) are used to compute C^CYjZ) and 

Now the computational plane and its boundaries are completely defined 
so that for any mesh point (X,Y,Z) in the computational plane the corresponding 
point (x,y,z) in the physical plane and all the necessary transformation de- 
rivatives can be computed. 


H^(X,Z). 
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z= f=Z 

Figure 6. - Three coordinate systems/ spaces usedo 


- 20 - 





Figure 7. - Body mesh point distribution. 





Figure 8. ■- Grid lines (near the body). 
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PHYSICAL MAPPED COMPUTATIONAL 

SPACE SPACE SPACE 


Figure 9 - - Plane Z = Constant in the physical, mapped and 
con^jutational spaces. 



Figure 10. - Series of mappings. 
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COMPUTATION OF REGIONS OF CONTINUOUS FLOW 
In the physical plane (xjy,z) the Euler equations are: 

wP + w = -(uP + vP + Yu + Yv ) 
z ' z X y X Y 

wu = -(uu + vu + TP ) 
z X y X 

wv = -(uv + w + TP ) 
z X y Y 

TP +WW =-(uw +vw) 
z z X Y 

wS = -(uS + vS ) 
z ' X Y 

where T = T/T^, P == ln(p/p ), S = (S-S )/c and all velocities are non- 

OO ^ ^ ® ^CO 

dimensionalized with respect to Poo/poo (the barred quantities are dimen- 
sional), X = x/l, Y = y/i-j z = z/t (Z is an arbitrary length), 

The equation of state for an ideal gas becomes : 

ln(T) = P (y-1) + S (20) 

Y Y 


(a) 

(b) 

(c) 

(d) 

(e) 


The dependent variables are P, S and the Cartesian velocity components 
u,v,w (Fig. 4). Now transforming all derivatives to the computational plane 
the following results are obtained. 


where f is the vector 


f 

= f^X + f^Y 

+ t^Z 

(a) 

X 

X X Y X 

Z X 

f 

y 

= f^ + f Y 
X y Y y 

+ f„Z 

z y 

(b) 

f 

= f^X + f^Y 

+ f z 

(c) 

z 

(P,U 

X z Y z 

,v,w,S) and 

Z z 


( 21 ) 
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(a) 


X - X r 
X r X 


+ x„e 

0 X 



X 

y 


X r + x^e + x„2 
r y 0 y p 


X 

z 


X r 
r z 


+ X e 

0 z 



(t) (22) 

(c) 


Similar expressions can be written for the Y and Z derivatives. The deriva- 
tives of (X, Y,Z) with respect to (r,0,2) and (r,0,2) with respect to (x,y,z) have 
already been discussed in the previous section. 

Combining equations (l9a) and (l9e) the following form of the Euler equa- 


tions is obtained which are used in the present solution. 

"z " -^"■2lV^22'b(^ ^21^Y^ ^22^Y^ 

■ (23) 

^Z "" ■^^31^x'^^33^x'^^31^y'^^33^y^ 

^z = -(^i''x^%2^x^%3^x^^44^x^Vy-'\2"y-^\3^y+^4^y) 


where the coefficients appearing in equations (23) are defined in Appendix B. 


At a data plane, Z = Z^ = constant, all the quantities on the right side of 
equations (23) are known and therefore the derivative f^ can be computed and 
used to predict the dependent variables at Z = Z^ + AZ. 

The step size AZ in the marching direction must satisfy the C.F.L. condi- 
tion for stability (ref, 13)- If \y, are the characteristic slopes in the 
X,Z plane and are the characteristic slopes in the Y,Z plane the 

stability criterion is written as follows: 


AZx_ = AX/a^„ 


AZy^ = AX/Xy 
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Each of these quantities is evaluated for every mesh point at the station 
Z = Zq, and aZ is taken as 70^ of the minimum of all of these aZ values. 

A modified MacCormack, two-level, predictor-corrector finite-difference 
scheme (ref. l4) is used to integrate equations (23). It can be proven that 
the MacCormack scheme is accurate to second order for a linear system of 
equations. So that the truncation error is of the fom. 


where t is a length in the physical plane. In regions where is large 

the mappings tend to assure 0 so that the truncation error remains small, 

while keeping the total number of grid points to a minimum. 

Equations (23) can be written in the following general form 

f„ = [A]f^ + [B]f,. (25) 


where as previously defined f is the vector (p,u,v,w,S) and [A] and [B] are 
matrices of the coefficients of equations (23). With these equations the 
MacCormack scheme proceeds as follows. 


Level one ; 


f„(Z ) = [A]fv + [B]fv (all quantities evaluated at Z ) 

Zf O A I . O 

= f'(Zo) + predicted value) 

Level two : 


(a) 

(b) ■ 
(26) 


'f = [A] f^ ^ [L]fv (all dependent variables evaluated with (c) 

the predicted values and all independent 

variables are evaluated at Z = Z + AZ) 

o 


f(Zo +Az) - (f + f(Z^) + ? 2 AZ )/2 


(d) 


The f^ derivatives are taken one sided in the positive X-direction in 
level one, the f^ derivatives in the positive Y-direction. For level two 
the direction of these derivatives is reversed. 
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This procedure defines all the dependent variables at interior points 
of the computational plane ( 1 <NN<WG(l) and 1<1V1M<MC ( I ) ) ; figure 11. The body- 
point calculation and the shock point calculation -will be discussed later. 
However, note here that all imbedded shock points have two mesh points 
associated with them, one for the low pressure side and one for the high 
pressure side, both having the same position in the physical plane (Fig. 11 ). 

The low pressure side of all shocks are computed following the MacGormack 
scheme and taking X and Y derivatives into the low -pressure region in both 
levels. The low pressure side of the bow shock ((M =NC(L), L = LC)) is- 
defined by the given free stream conditions. ' - 

The points on the symmetry planes (MM =1, I == 1 and MI'l - MC(l), I = IC) 
are computed using the same scheme and the symmetry conditions P.^ = v.^ = w^ = 

Sy = 0 and u = 0. The points on the internal boundaries that are not shock 
poin-bs are also computed using the MacGormack scheme. In level one the points 
KN=1 and MM=1 are computed, taking the difference between M=l, NN=2 and MM=1, 
MM=2 for the X and Y derivatives respectively. After level one quantities at 
NG(l) and MC(L) are updated (i.e., f(WG(L),M). = f(l,M),. ,, and f(N,MG(l)),. = 
f(W,l)j_^^). In level two the -points on the other side of the surfaces M = 

NG(L) and MM = MG(l) are computed and afterward the points M=1 and MM=1 are 
updated. 

The modifications to MacGormack' s integration scheme were originally found 
necessary in the calculation of blunt nose bodies. For blunt nose bodies 
derivati-ves aci-oss the "edge" of the entropy layer are not allowed. The entropy 
layer calculation will be discussed later, for now assume the position of a 
surface X = F(Y,Z) representing, in the computational plane, the edge of the 
entropy layer, (Fig. I2a) and all dependent variables on this surface are known. 
Across this surface the derivatives S^, u^, v^^ and w^ become very large as the 
surface approaches the body. In computing these derivatives at the mesh point 
N of figure 12b, instead of taking the X derivatives between points N and N + 1 

differences are taken between W and the * point which is on the entropy layer 
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surface. If the entropy layer surface point becorn.es very close to the mesh 
point (Fig. 12c) the dependent variables at N are set equal to those computed 
at the entropy layer surface point since the distance between N and * is too 

small to compute S„, u„, v^ and 'w between them. 

A X A X 

In order to calculate the flow field over bodies with blunt nose entropy 
layers it was found necessai’y to use "windward" differencing. Equations (23b, 
c, and e) state the variation of u, v and S along a streamline^ i.e., the 
velocity direction is the characteristic direction of these equations. Accord- 
ingly, for the derivatives v^, u„ and S in these equations, differences are taken 
in a direction determined by the velocity direction in both levels of the 
MacCormack scheme (Fig. 13) J the same is true for v^, u^ and S^. If the slope 
and/or "the velocity vector is small the derivative direction is changed 

between level one and two as usual. Since in equations (23b, c, and e) infor- 
mation is carried along streamlines, windward differences satisfy the rule of the 
domain of influence, so that windward differences are used even when there is no 
blunt nose entropy layer in the computation. The techniques of following the 
entropy layer and using windward differences were originally suggested by Moretti 
and Pandolfi (refs. 7 and 8). When derivatives are approximated with windward 
differencing they are no longer formally second order accurate. But it was 
found that windward differencing yields more stable results and the integration 
scheme was found to be second order accurate in a numerical experiment. 

The calculation of interior points is the most time consuming part of this 
computation mainly because it is done many times. The scheme used here keeps 
the computational time to a minimum by keeping the total number of mesh points 
as small as possible. 
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MM 
M = 



Figure 11. - Region and mesh point notation. 



Figure 12. - Entropy layer surface (computational space 



Figure I 3 . - Wind^vard differencing scheme. 




TREATMENT OF SHOCKS 


In this section the computation of grid points on the high pressure side 
of all shock waves, the detection of imbedded shocks and the intersection of 
two same family shocks will be discussed. In this area we draw systematically 
from Moretti's extensive research on the treatment of shockwaves (ref. 1 and 
6 ). 


The bow shock, all wing type shocks (i.e., imbedded shocks which in 
general originate near the body and move toward the bow shock, caused by 
canopies, wings and vertical tails) and cross flow shocks in the flow field 
are computed as discontinuities satisfying the Rankine-Hugoniot conditions. 
The bow shock and wing type shocks are defined in the mapped plane by 
r = c(9,p) and the crossflow shocks by 0 = h(r,^). At a data plane all 
dependent variables are known, and in addition, the quantities c, c^, c^. 


h, h^, h^ are also known. Using equations (l6)to(l8) C(Y,Zq), Cy,C 2 ,H(X,Zq), 
and are confuted for all the shocks in the flow field. At Z^ + AZ the 


positions of the shock points can be computed by using; 


C(Y,Zq + aZ) - C(Y,Zq) + G 2 (Y,Zq) aZ (a) 

(27) 

H(X,Zq + aZ) = H(X,Z^) + Z^) aZ (b) 

Then Cy(Y,Z^ + aZ) and Hj^(X,Z^ + AZ) are computed using central differences. 

With these quantities and equations (l6) c, c , h and h can be congjuted at 

9 ^ 

Z = Zq + az* 

After the first level of the MacCormack scheme the predicted values 
of the dependent variables on the low pressure side of all shocks are computed 
(the variables on the low pressure side of the bow shock being the constant 
free stream values). With the low pressure side of the shocks known the high 
pressure side is computed by an iterative process. 
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A value of h or c is guessed, between the values corresponding to an in- 

S" ir 

finitely weak shock and the value which gives a subsonic axial Mach number. Once 
this guess is made the normal to the shock can be defined. Let 


F s r - cCe,^) 

or F = e - h(r,^) 


Then the normal to the shock I is given by: 


' / ^ O O ^ ^ 

I = (Fi + Fj + F k)//F +F +F =I, i + I..,j + I.,k 
yi .Y z X y z 12 3 


where 


F = Fr + F0 + F 2 
X r X 9 X 


F = F r + F^e + F^2 
y r y 9 y 


F = F r + F^0 + F^2 

z r z 9 z P^z 


With the normal to the shock defined, the Rankine-Hugoniot conditions can be 
applied. Using the subscripts 1 for the low pressure side and 2 for the high 
pressure side we have; 


V = V • I 
^ 1 


M = V /vA^ 
n^ n^' ' 1 


pg/pi = (y+1)/2.]/[1 + (y-1)/2.] 


(a) 

(b) 

(c) 


= [(y+1)/(y"1)p2/Pi -1-1/[(y'^iV(y-1) - P2/pi3 

= [Ti(P2/Pi)]/(p2/Pi) 


% = V - V I 
2 1 n^ 


(e) 

(f) 


( 28 ) 
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Where M and V are the Mach number and velocity normal to the shock and V 
n n J- 

is the velocity tangent to the shock. 

An intrinsic coordinate system is defined at the shock with the three 
directions (IjJjK), coordinates (§,T),(b) and velocities (iTjVjw) such that I 
is normal to the shock and: 

K = (i X k)/]! X kj = K^i + + K^k 

J = i X K = i + 1^3 + J3k 

In the HjU) plane the characteristic that intersects the shock from the high 
pressure side has a slope 

^_^_(uiT+a /u"^ + w ^ - a^ ) (29) 

' ' (»2 - ab “ 

Figure l4 shows the shock, characteristic (in the plane) and the point 
(*) in the Zq data plane where the characteristic originates. The characteristic 
slope at the shock point is first evaluated and then the position of the (*) 
point is computed using the relations: 

(„* = - aZ/(K2 + U3) 

I* = 1 cu* 

y"" " ^SH ^ 

Where the subscript SH refers to quantities at the shock at Z + (Fig. l4). 
Dependent variables at the (*) point are obtained by linear interpolation. 

A value of the pressure on the high pressure side of the shock is computed 
using the compatibility equation along the characteristic: 
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~ 2 /p /~ 2 ^ ~2 2 / 

'?, = yw/'ia./vL+v - a. (aj 

\ - (u w + a + w ^ - a^)/ (w ^ -a^) (b) 

R = [(u - Xw) (v + yv^) - Y V ) Y V w^1 (c) 

^ /~2 ^ ^2 2 

a / u + w - a / . 

( 31 ) 

dT = uV w* - Wgjj (d) 

= P* + Ru)* - gdT (e) 


Where p,!, and R are averaged between the * point and the shock point. 

The iteration is continued until this value of pressure agrees with that 
computed from the Rankine-Higoniot conditions, for some value of c^ or h^. 

For weak shocks this iteration may converge to a pressure ratio P 2 /Pj^ <lj 
in these cases the value of c or h which gives Pp/p- = 1 is taken. 

Cross flow shock points at the body must satisfy the body boundary condi- 
tion, that is the velocity normal to the body on the high pressure side of the 
shock must vanish. This implies that the shock normal, at the body, must be 
perpendicular to the body normal. This condition gives a relationship between 
h and h at the body: 

h = (F„ - h F^ )/F ( 32 ) 

r H p Bz ' r ' 

Where : 


F„ = F„ 9 + F„ 9 + F 0 

H Bx X By y Bz z 



F r + F r 
By y Bz z 
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and , F_ are x, y, z derivatives 

Bx By’ Bz 

of F„ = r-b(n,y) 

D 

and r = b(9,2) defines the body. 

After the second level of the MacCormack scheme the corrected, final 
values of the dependent variables on the low pressure side of the imbedded 
shocks, the values of c and h„ confuted after the first level, and the Rankine- 
Hugoniot conditions are used to compute the final values of the dependent 
variables on the high pressure side of the shocks ^ 

The first problem encountered when one treats imbedded shocks as dis- 
continuities is their detection. There have been a number of techniques pro- 
posed (see ref. 6). One of the earliest procedures has been found to be well 
suited for the type of shocks encountered in this problem. 

Cross flow shocks and wing t3^e shocks are detected in very similar ways. 
For cross flow shocks, the pressure distribution is monitored in the Y- 
direction and for wing type shocks the pressure distribution is monitored 
in the X-direction. At a data plane Z = the maximum pressure gradient 
for wing type shocks and for crossflow shocks is located. Then a third 
order polynomial is fit through the four mesh points adjacent to the maximum 
gradient (Fig. 15)- This polynomial takes the form: 

X = a^ p3 + a^ P^ + a^ P + a^ (33) 

Where (X = X for wing type shocks) 

(X = Y for cross flow shocks) 

and the coefficients a^, a,^ and a^ are computed by matching the curve fit 
to the four mesh points (Fig. 15). The condition used to determine the origin 
of a shock is dx/dP = 0 vmich implies dP/dX “* <». Applying this condition to 
equation (33) yields an equation for P^ of the form: 

’’f ° -«i 1 A/ - 3 Va 

a 

o 
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When 




- = 0. This equation has one real root. When this condition 


3 2 

is satisfied a shock is inserted in the flow field at X„ == a P-:. + a, P,, + 

f of 11 


‘‘ah ^ 


3’ 


Cross flow shocks are assumed to originate on the body, so that the 
pressure distribution is monitored in the Y-direction on X = 0 (the body). 
Once a shock is found on the body monitoring is begun at increasing values 
of X = constant. 


In general it is not known at what value of Y the first shock point on a 
wing type shock will be found so that the maximum pressure gradient P^ at all 
values of Y must be tested until the first shock point is detected. Once a 
wing type shock is detected additional shock points are sought at the grid 
points adjacent to the end shock points (Fig, l6). 

Finally, consider the intersection of two shocks of the same family'. Cross 

flow shocks do not move very much (i.e., h. is small) so that a scheme to handle 

r 

the intersection of two crossflow shocks is not needed. Wing type shocks are 

detected near the body (they are usually caused by compressions on the body) 

and they move toward the bow shock, so that they are all of the "same family". 

As a shock moves toward the bow shock the region between the imbedded shock and 

the body gets larger and the region between the imbedded shock and the bow shock 

gets smaller so that mesh points must be moved from the outside region to the 

inside region. When the distance between two shocks becomes a small percent 

of the total shock layer (l to 5^) at some value of Y (Fig, 17) the two shocks 

are intersected. A local, exact, two dimensional calculation is used to conpute 

the intersection of two shocks. The same intrinsic cqprdinate system is used as 

was discussed in the shock point conputation, where I is normal to the outside 

shock. In the §,«) plane the intersection of two same family shocks is shown 

in Fig. l8. The conditions in regions 1 and 4 (Fig. l8) are known. The slope 

c of the resulting shock is assumed and Rankine-Ifiigoniot conditions (28) and 
r 

the conditions in region 1 are used to compute the conditions in region 2. The 
pressure in region 3 is set equal to that in region 2 (since the pressure is 
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constant across a contact surface) and the total pressure in region 3 equal 
to that of region 4 (since the total pressure is constant through an expansion 
fan). With the pressure and the total pressure in region 3 the Mach number 
can be computed. Then from the Prandtl-Meyer expansion relation the flow 
direction 6^ is obtained. The iteration is continued until the velocity 
direction 6^ matches that computed in region 2. 

All the iterations follow the same procedure. Assume we have two 
functions of a single variable g(Tl) and G(ti). The problem is to find the 
value of p for which g = G. Two values of T] are assumed and two errors 

~ ^^2 ~ ^2 ~ ^ 2 ~^ computed. With e2L5e2>''l3_ \ ^ linear 

variation of e vs T] is assumed to predict the value of p which will force 
e — ^ 0 • 

Tl^ “ ilp " (II 2 - n^)/(e2 “ Gp) 

The assumed value is repeatedly updated in this manner until g "* This 
scheme was found simple and fast and in most cases it converges in 4-5 
iterations. . • • 


- 35 - 



X 


^SH' ^SH 



Figure l4. - Shock point calculation. 



Figure 15 . - Four mesh points used in shock detection. 


ADJACENT VALUES 


END SHOCK 
POINTS 


Figure 16. - End shock points (wing type shocks). 
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17. - Shook Ihtersectlon 1 , 


constant plane. 
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SHARP LEADING EDGE SHOCKS 


When a confi^ration has a sharp leading edge wing (Fig, 19) its wing 
leading edge is ccanputed using a local two dimensional solution. The component 
of the Mach number normal to the leading edge must be supersonic in order for 
this technique to be applied. Also the wing must be sharp from the root. 

The mesh point distribution near a sharp tip is shown in figure 20, 

The mappings are not modified for this type of configuration but a cross flow 
type surface is inserted starting at the wing tip, in order to have a double 
point at the tip, on the top and bottom of the wing. 

The shock points at the tip (the case of an expansion fan on one side 
and a shock on the other can be handled) are automatically inserted when the 
wing starts. All other points on this shock are detected as in the case of a 
blunt nose wing. The extensions of the shock are treated as discussed previously. 
If the flow direction is such that there is a centered expansion fan on the top 
or bottom of the wing (Fig, 20), the expansion is computed explicitly at the 
tip point, and an arbitrary surface is used from the tip to the symmetry plane. 
Thus, the expansion fan is computed using the finite difference scheme for all 
points except the one at the tip. 

The calculation of the tip points utilizes an intrinsic frame of reference 

( §,T|,a)), (u, v,w) defined in figure 21. The T1 direction is tangent to the sharp 
edge of the wing (hence, also tangent to the shock) so that 7 remains unchanged 
across the shock. In the §,u) plane, a local two dimensional wedge calculation 
is performed. The conditions on the low pressure side of the shock are known. 

A third order algebraic equation (ref. 15) is used tcJ compute P, the shock 
angle, and then the Rankine-Hugoniot jump conditions are used to compute the 
dependent variables on the high pressure side of the shock. If the flow is 
such that there is a centered expansion fan at the tip point the Prandtl- 
Meyer expansion fan equations are used to compute the dependent variables on 
the body surface. This calculation is done after each level of the McCormack 
scheme . 
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This method will not handle the situatioxi in which the lVfe,ch number normal 
to the leading edge of the wing is subsonic or the situation in which the 
wing angle ( 5 , Fig. 21 ) is large enough to force the shock to be detached. 
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Figure 19. - Sharp leading edge wing configuration. 


SHOCK 


(a) (b) 

Figure 20. - Mesh point distribution, sharp leading edge wing. 



|3 


Figure 21. - Sharp leading edge wing intrinsic frame. 




_4o- 




BODY POIM COMPUTATION 


The boundary condition at the vehicle surface is u = 0, where u is the 

velocity normal to the body. The entropy at the body is computed by using 

equation (23e), as for any other mesh point. At the body the coefficient 

of S (i.e.j a ) in equation (23e) vanishes, so that this derivative does 
X IPP 

not affect the calculation of entropy on the body. 


To compute the pressure on the body the continuity and three momentum 
equations (23a, b, c, and d) and the body boundary condition are com- 
bined to write a compatability equation along the characteristic (in the X,Z 
plane) reaching the wall from the flow field. 


^ ^'ii^x ^ ^12 ’^x ^13’^x ^l4 ^X ^1 


(a) 


Ur 


‘2i.^X ^22 ^X ^2 


'^Z ®'31^X ^33 ^3 


wz + + a^^2 '^X^ ^43 ^ ^44 "'x = ^4 


(b) 

(c) 

(d) 


(34) 


where 


^1 = -^\l^Y ^ ^2"y ^ \3^Y ^ \4^Y^ 
^2 "" " '^21^Y ^22^^Y^ 

R3 = - (b^^Py + b^3VY) 


1^4 " - ^^4i^Y ^ ^ ^3^Y ^ V^y) 

Substituting for a^^? ^ equation (34a) becomes 

^z 1 ^x '^x ^1 


(35) 
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where : 


T = u/w and a - v/'w 

Taking the difference between the product of w with equation (3^b) and u 
with equation (3tb) and substituting for ^43’ ^42 ^44 following 


result is obtained 
2 


w + P^^(wa23_ - uaj^^) + w^(X^ + X^tw^ - X^o) 


'4l^ 


X 


(36) 


+ w^ryT = wR^ - 


Similarly, taking the difference between the product of w with equation (34c] 
and V with equation (34d) and substituting for ^■43’ ^-42 ^44 

following result is obtained 


^"®31 ‘ '^41^ * 3l!L * V^“x 


2 ^1 
^ O vT X^Tj^ = wR^ - vRj^ 


(37) 


The body boundary condition is u = 0. Since X = constant is the body, this 
boundary condition can be written as 


u = uX + vX + wX =0 
X y z 


or 


tX + aX + X =0 
X y z 


(38) 


Combining equations (35), (36), and (37) and using equation (38) the 
equation for the characteristic slopes can be written in the form 


^ 1 - -a. 4T + / a ^t 2 + ^ 

A. 

w w^ 1 


W 


(39) 


X is the slope of the characteristic reaching the w'all from the flow field 
and the compatibility condition along this characteristic is: 
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( 40 ) 


[X - yT (uX.^ + vXy)] (P^ + XF^) 


Kv 


* »!y (X^Tj. + XyO^) + (X^T^ X^a^) = E 

Ap ^ 




where 


R - /x - YT (uX + vX ) \ R 


ApW 


X y 1 


+ X,^Y^ ^wR 2 - uR^) + X vX (wR^ - uR|^) 


The equation for the body can be written in the form 

F = r - BCi^Z) 

Thus the body boundary condition is 

xF + cjF + F =0 
X ^ y z 

This equation holds for all values of Z, thus 

a J + T J = - (F _ + aF „ + TF „) = R 
ZTy ZT X zZ yZ xZ 

Wow using equation (lOa) equation (4l) becomes 

F = r - B = (C^ - B) X 

Thus, differenting equation (HU) 

F = [(C.) - B ]X + (C. - B)X 

X t X X-" t X 

F = [(C) - B ]X + (C, - B)X 

y ^ I y y-‘ -I ' y 

and for X = 0 


F 

X 

F 

y 


("t - ®)^x 
(C^ - B)Xy 


(Ul) 

(U2) 


(U3) 

(UU) 


(a) 

(b) 


(U5) 


(a) 

(b) 


(U6) 
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Using equation (43) and (46) in equation (4o), the following result is 
obtained 


P = 
z 


R Yv£ 
1 + ■ 


_J 


+ 


(X 
X X 


X a ) 

y X _ 




(4t) 


[A 


+ 




This equation is integrated with the same scheme used for interior points with 
the X-derivatives computed using three-point end differencing away from the 
body. 


To compute the velocity components on the body an intrinsic frame 
(I,J,K) is used with velocity components (tT,’v,w). The vector I is the unit 

A A 

normal to the body with J and K defined as follows: 


/V ^ 

I = II i + + 13 k 


J 


(I X k)/ 


I X k = 


i + JpJ 


H- 13k 


(a) 

(b) (48) 


K=ixK=K^i+K2j+K3k 


(c) 


A vs . 

where i, j and k are defined in figure 4. The x and y momentum equations are 
used to compute v as follows. Equations (34b) and 34c) are integrated using 
the MacCormack scheme to obtain u and v and then v is obtained from the 
equation 


V = uJ^ + vJ^ (49) 

From the integrated form of the energy equation the ^ component of velocity 

can be obtained . — 

w = / 2H - 2YT - v^ ( 50 ) 

O r 

Y-1 

where T is computed from P and S and the equation of state (l9)- The three 
Cartesian velocity components are 


u = 




(a) 


V = V J2 + w 
w = V J 3 + w 


(b) (51) 

(c) 
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Thus, all the dependent variables (P,u,v,w, and S) on the body are defined. 

Modifications of this calculation for real gas and entropy layer 
effects will be discussed in later sections of this report. 



BLUM NOSE EMROPY LAYER CALCUIATIOW 


On blunt nose bodies, as the computation proceeds downstream, the entropy 
gradient normal to the body becomes very large. Stream lines that cross a 
weaker bow shock and therefore have low entropies approach the body which is 
wetted by the stagnation stream line and therefore has a very large entropy. 
The pressure gradients at the body remain small as the edge of this layer 
approaches the body but the normal derivatives of S,u,v and w become very 
large. This physical phenomenon can create numerical problems for a calcu- 
lation which does not handle it properly. A technique proposed by Moretti 
and Pandolfi (refs. 7 and 8) is used to account for this phenomenon. 

In this section, after defining the "edge" of the entropy layer, the 
detection of points' on the edge of the entropy layer, the calculation of 
the dependent variables at these points and the modification of the body 
calculation when the edge becomes close to the body will be discussed. 

In figure 22 the entropy distribution (S vs X) in the windward symmetry 
planes on two typical geometries at several values of Z (axial stations) 
is shown. The *'ed points in figure 22 denote what is called the edge of the 
entropy layer. The entropy distributions are similar in other circumferential 
planes up to the top symmetry plane. The locus of these * points defines an 
"entropy layer surface" as r = R,^(Y,Z) (Fig. 23). This surface originates 

nij 

at the bow shock and moves toward the body as one proceeds downstream. 

This surface is not a coordinate surface, so that in the computational 
plane the * point will be between two mesh points (Fig. 12) at each value 
of Y (it should be noted that this surface does not originate at the same 
axial station for all circumferential planes). 

The key idea of this procedure is that no derivatives should be taken 
across this surface. The method used to insure that the derivatives S u„, v„ 

and w^ are not taken across this surface when computing the mesh point near 
it is discussed in the section on interior point computations. The questions 
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remaining are : how the surface is detected and traced and how the dependent 

variables are computed on it. 


The surface is detected at it's origin, the bow shock. When the entropy 
distribution is similar to the one shown in figure 22a (the entropy has a 
minimum) the * points at each value of Y is located at the mesh points adja- 
cent to the bow shock when the entropy there is a minimum. The surface at 
this Y is initiated with the value of the dependent variables at this mesh 
point and is tracked separately from this station on. When the entropy dis- 
tribution is similar to the one shown in figure 22b (no minimum exists) the 
* point at a value of Y is started when the S derivative at the mesh point 
adjacent to the shock has a minimum (i.e., S sO). 

K/j 

The surface r = R^(Y,Z) is a stream surface (i.e., a surface containing 

xlL 

the same group of streamlines for all axial stations). This means that the 
velocity normal to the surface is zero. When at a station Z, all the dependent 
variables on the surface and the position of all the * points are known. 
Therefore, '^^e surface can be computed using finite differences. 

The normal to the surface is 




( 52 ) 


Since the velocity normal to the surface vanishes 
From this equation the following result can be obtained 


X x ' Y HLy' x^ 


The position of the entropy layer surface at Z^ + AZ is given by 




( 53 ) 


( 54 ) 


( 55 ) 
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Since the pressure is continuous across the surface the pressure at each 
* point can be obtained by interpolation using the two adjacent mesh points. 
The entropy and crossflow velocity at each * point are computed using the 
following relations 

= - [Sy, (uY^ + vYy + wY^)]/w (a) 

U2 = - [T(Py^^ + ^ (56) 

V2 - - [T(Py\ + P X^y) ^ '^^y ^ 

where S ^' 3 and v^’ are the Y-derivatives on the entropy layer surface. 


These three equations are derived by taking a coordinate system (X',Y',Z') 
defined by the following transformation 


X' - (r - B)/(Rj^ - B) 

(a) 


Y' == Y 

(b) 

(57) 

Z' = Z 

(c) 



The substantial derivative of any quantity f = (u,v, S) in this coordinate 
system is 


Df 

Dt 


= uf + vf + wf 
X y z 

= f, ,(uX' + vX' + wX') 

X X y z 

+ f^, (uY' + vY' + wY' ) 
Y ' X y z 

+ f^, (uZ' + vZ' + wZ' ) 
Z ^ X y z 


(58) 


The coefficient of zero on the "entropy layer" surface since the 

velocity normal to that surface is zero, while f^' is the Y-derivative on the 
surface and 


■ru. ^ 

JJU 


V' 


(uY + vY 


wY ) 
z 


fzW 


(59) 
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To compute the velocity components an intrinsic coordinate system is 
used with unit vectors (I,J,K) defined as follov^s: 


i = (x i + X j + X k)/ /x^ + X^ + X^ 
X z " ^ y. y z 

(a) 


J = (I X k)/ 1 i X k 1 

(b) 

( 60 ) 

K = 1 X J 

(c) 



A. ^ ^ ^ ^ 

The velocity components in the I,J,K directions are called iTjVjW. The 
streamline slope in the (X,Z) plane is continuous across the entropy 
layer surface (just as in the case of a contact surface, where the slope 
of the contact surface is the same as that of the streamline adjacent to it 
ref. 7 and 8). Therefore (iT/w) can be interpolated from the adjacent mesh 
points. Moreover, v = uJ^ + vJ^ where u and v are computed from equations 
34b and 34c. Next, w and u are computed- from the following equations 

w = - 2YT/(Y-1) - V ^ (6l) 

(l+(u/w)2) 

u = (u/^Ow ( 62 ) 

^ ^ ^ ^ 

Then, from u,v,w and the vectors I,J,K,the velocity components u,v, and w 
are computed. 

When the entropy layer surface becomes very close to the body (that is 
within 1% of the total shock layer thickness, figure 23 c and 23 d) at some 
value of the circ\unferential coordinate the entropy and the crossflow velocity 
Cv) on the body are changed to their values on the edge of the entropy layer. 
The pressure at the body remains unchanged. On the body there is a jump in 
entropy at the points v/here the entropy layer surface separates from the 
body (Fig. 23c), so that Y-derivatives ca,nnot be taken across these points. 
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For grid points that have collapsed Y-derivatives are always taken onto the 
entropy layer surface. For the Y-derivative at a grid point that hasn't been 
collapsed but is adjacent to a point that has been, a pseudo grid point is 
located at the adjacent collapsed point with the blunt nose entropy, and 
velocity direction and pressure of the collapsed point. 

Changing the body entropy, from the high stagnation-point value to 
the low value at the edge of the entropy layer increases the axial Mach • 
number and therefore increases the step size (aZ) at which the ' computation 
can proceed. In seme cases the axial Mach number would become subsonic with 
the stagnation entropy and the computation would be unable to proceed. This 
significant advantage could not have been gained without following the entropy 
layer explicitly. 
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REAL GAS EFFECTS (EQUIUBRIUM/FROZEU AIR) 


In some flight regimes encountered by hypersonic vehicles the simple 
assumption of air being a perfect gas is not acceptable and real gas effects 
must be taken into account. Grossly speaking, real gas effects shift but do 
not change the qualitative pattern calculated with perfect gas thermodynamics. 

A code for real gas should, therefore, be patterned along the same lines as a 
code for perfect gas. 

It is in this spirit that the real air problem has been approached. 

The salient points of the procedure that has been adopted are: 

1. The equations of motion are written using the same dependent 
variables as for the perfect gas calculations, namely, the three 
velocity components, the logarithm of pressure and entropy. 

2. Simple analytical expressions, developed by Moretti (ref. l6), are 
used to fit the Mollier diagrams for equilibrium air. The ability 
to suddenly "freeze" the fluid with an equivalent ideal gas y 

is included. 

3. The computational time has not increased out of proportion. The 
computational time needed for a perfect gas computation is increased 
by a factor of less than 50 ^ when equilibrium real gas effects are 
included. 

4. An accuracy of around 5^ is maintained for enthalpies less than 
500 X (psi/PgL* ■ 

5. The parts of the code dealing with gas properties are kept separated 
from the flo\/ calculations so that if necessary the code may be easily 
extended to helium, nitrogen, CFj^, or supersteam. 

The equations of motion in nondimensional form are: 
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(a) 


wP + uP + vP + r(w + u + V ) = 0 
z X y z X y' 


ww +UW + w +tP = 0 
z X y z 


WU +UU +VU +Tp =0 

z X y X 


wv +UV + w +TP = 0 
z X y y 

wS +uS +vS =0 
z X y 

2 

where t = p/p and F = a /t. The only difference from perfect gas is in 
r and T. The quantity, t is not the temperature except for a perfect gas; it 
is T = zT, where T is the temperature and z is the compression coefficient. 
Similarly, F is not Y(cp^ /^^ao )®^^ept for a perfect gas. For equilibrium 
flow both T and F are obtained at every point via the Mollier fits. 

The Mollier fits for t and F developed by Moretti (ref. l6) are written 
in terms of enthalpy and pressure. 


(b) 

(c) 

(d) 

(e) 


and 



r = r„ = fe 
r = 


+ bh + c/(h + d) 

(0 < h ■< 50, 150 < h < 350) 

(a) 


+ PggOi-hn)^ 

(50 < h < 150) 

(b) 

i6h) 

+ ke^^/50-10) 

(350 < h < 500) 

(c) 


+ § h + Tle^^ )/H^ 

/ 0 < ii < 50,\ 

(l50 < h < 500 ) 

(a) 


+ Pgg(h + m)^) 

(50 < h < 150) 

(b) 

(65) 


where 
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and 


t = 12030.872 -.?64759(P + 157.7555) 

77.938126-P 

I = 12.813 + .0871477 (P + 9.4423)^ 
d - 250 ii - t)/(2t - I - 79.4) 

a = (.3176 - b) d 

b = 3176 + (1 + d/250) (r - 79.4)/250 

c = -ad 

f=-1.83^1.098/(e--^^S6 ^^-l- 72725),,_) 

g = -.0038 - .00219476/(P- + 10.3635) 
m = -84.'6 - .34744 P (1.- .21T45P) 
k = 9.2217 - .052l317l(P + 8.0605)^ 
e = 1.0459 + .00424528P 
§ = (12.707 - .424528 P) 10“^ 

T1 = 1.1828-e 

a = -(.001955 + 1 )/T\ 

The subscript SL refers to sea level conditions as follows (the bar in the 
following equations means dimensional quantities): 


^SL 

PSL " PsL^Pco 


( 66 ) 


The normalized entropy, S, as a function of enthalpy and pressure is given by 


S = (y --1) S/(p^ 


(67) 
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where 



4.82068 In (h) + 11.875 + .0245 h + 

R 

1175. /(h + 50) + .434294 B P 

and 

B = - 2307.- (.0042 h - .092)/[l. + ~ 

Since and p must be obtained as functions of entropy and pressure, the 
expression given for the entropy is solved for the enthalpy, and 
then with the enthalpy and pressure, t and F are obtained from the 
remaining fits. 

The interior point computation can be performed by the same procedure 
used for perfect gases, except that t and F are computed via the curve fits. 

Some changes are necessary in the shock computation. If subscripts 1 

and 2 refer to the low pressure side and high pressure side of the shock 

— ♦ 

respectively, and if and refer to the velocity components in the 
directions of the normal and tangent to the shock respectively then the 
jump conditions are given by 



(a) 

V 1. + 

(b) 

Pl/p2 = ^n2/^nl 

(c) 

P, /P^ = ^ ^ ^2^2 

(d) 

1 + f^M^ 

(e) 

^2 ^ 2 



( 68 ) 


*Note : 


S is the dimensional entropy referenced to 0° Kelvin and 1 atmosphere 
of pressure. 
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The method of solution follows the same procedure used for a perfect gas. 

The shock slope and V are guessed and equation (68a) is used to solve for 

^2 

^2^^^ . Equation (68d) is then solved for and from equation (68e) hg 
is obtained. The ratio pg/p^ is then computed from equation (68c). 

A second value for pg/p^ is obtained from the curve fits of the Mollier 
chart using h^ (obtained from Eq. (68e) and Pg. An iteration is now per- 
formed on until the difference in the two computed values of pg/p^^ is 
sufficiently small. The rest of the calculation is identical to that used 
for a perfect gas. 

The equilibrium air calculation starts at the nose of the vehicle and 

continues until a "freezing plane" (i.e., a plane Z = constant), after which 

both the chemistry and vibrational relaxation are assumed fixed (this is called 

frozen air in what follows). This freezing station (Z. ) is user input. At 

fr 

the freezing station all the dependent variables are related by the Mollier 
fits and are in equilibrium. (P^^, S^.^, u^^, v^.^, w^^, and h^^. The 
frozen equations of state are (ref. 1?) 

Pfr/Pfr = ^fr (a) 



(b) (69) 

(c) 


with the quantities (R/K^), T and defining the properties of the gas 

downstream of freezing station. In order to match the equilibrium state 
(i.e., u^^_ = UgQ etc.) at each mesh point one must set: 
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£/R^ (Peq/^EQ^/'^EQ 

"eQ «/«„-’ic!> 

"^EEF - T^) ’kq " (r-1) 

Downstream of this station, these quantities are constant along streamlines 
because they are functions only of the concentration of species which are 
in turn constant along streamlines. But since the variation F and 

ScTTin is small at the freezing plane, they are averaged. 

AHji? 

Downstream of the freezing station, equations (70) are used instead of 
the Mollier-curve fits to relate thermodynamic properties. 


(a) 

(b) (70) 

(c) 
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SPECIALIZED OUTPUT 


Calculations have been added to the procedure which use the inviscid 
flow field solution to compute desired results in a useable engineering form 
or as an input to other calculations. At the user's option, the flow field 
data may be used to calculate aerodynamic coefficients, boundary layer input 
quantities, (streamlines, metric coefficients, etc.) and sonic boom data. 
Following is a description of the approach used for each of these calcula- 
tions. 

Aerodynamic Coefficients . Aerodynamic coefficients are calculated, when de- 
sired, along with the flow field. Vectored pressures, located at the area 
centroid of triangular facets (two such facets between each pair of mesh points 
in the previous and current data planes) and directed normal to the facet plane, 
are summed at each step over user- specified portions of the vehicle. By mak- 
ing use of the body lines (running cross sectional control points) in the 
QUICK geometry package, the coefficients C^, C^^, C^^, and C^ can then be 
computed from the appropriate direction components of the integrated forces 
and moments for virtually any body components (even multiply disjoint sections 
such as the body-alone represented by the solid line in Fig. 2k). 

The pitching moments are referenced to a user- specified position in the 
vehicle symmetry plane. Coefficients are expressed in two fashions - (l) re- 
ferenced to the surface area integrated thus far for the given component (this 
surface area is included in the output), and (2) referenced to a user input 
area. In both cases the reference pressure and density are p and pm. 

Aerodynamic coefficients and wetted areas are available in terms of the 
entire vehicle, its components, and their distributions in z. 

Boundary Layer Input Quantities . The boundary layer input quantities are calcu- 
lated in two steps. During the flow field computation, integration techniques 
are applied to obtain the metric coefficient h^, defined herein. 

Consider Figure 25, in which the fixed body coordinates (x, y, z), the unit 
vectors (i, j, k). and velocity directions (u, v, w) have been indicated. 

Also shown is an inviscid streamline on the surface of the body, and the asso- 

AAA 

ciated coordinates (^, Tl, Q) , unit vectors (§, T\, C)> and the velocity direc- 
tions (u, V, w) . The 11- direction is taken along the streamline, the 
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^-direction is perpendicular to the hody and n and the ^ - direction is perpendi- 
cular to both T] and 

Associated with these two sets of coordinates, (x, y, z) and (§, T], ^), 
we have the three scalar factors (h^, h^, h^), defined by: 

^2 2 ^ 2 ^ 2 

hi = \ 

^2 2 2 2 

^ * >'n * "t| 

1,2 2 ^ 2 ^ 2 

*■3 “ "C " ’'c " "c 


Wow. 


%= h^Tl = h2(u/qi + v/qj + w/qk) 


or. 


^T| = hg u/q 
Yti = ^2 v/<l 
= ^2 w/q 

The body surface noimal {Q = C^^i + £2^ C3^) readily available, and 


= ^3 C2 


A. A, 


and since § = T] x 


‘'b " 1 ("C2-vC3)/q 

y^ = h^(x,^ z^-z^ x^)/h^ ^3 " ^ (’^3-wCi)/<l 

^3 = 
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In three-dimensions, "between any two points P and P' , we have: 


dP = d5 + dTl + dC 

= h^ dd + hg dTlIl + h^ dCC 


We get : 




r) P 


K ^3 ^ 


and. 


STl d? ^ d5 


thus, 


d /dP\ ^ /dP \ 

dH ^dC^ ^d1l^ 


0 |> dP \ _o_ / or s 

d| ^dC^ " dC ^d?^ 


d?- 


dll ^ dT] 


dhp „ 

— h -n + ^ h 

d§ '' d| 2 


— - C + ^ K = — - 

dll ^ dll 3 dC 


11 


"I? "2 


dh„ 


c + 


K 


h^ = 


dh-, 

= I + ^ 


h. 


a§ ' d§ ”3 dC dC 1 

Finally, dot multiplying by the appropriate unit vector (5, T] or Q) we obtain. 
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(a) 


Sh. 


1 

BTl 


^ a? 




h ^ 
1 aT] 


(b) 


( 71 ) 


ah. 


at] 


3-iul? 


^ ac 


(c ) 


We make the usual assumption of = 1, which is valid in the absence of 
a strong entropy layer. Observing equation (71a) , one will note that h^ is 
differentiated with respect to (§) which leads to a multiplier of the h^ de- 
rivative that corresponds to a cross flow velocity which may vanish, thus mak- 
ing an axial marching technique inapplicable to the determination of . The 
other parts of equation (71) involve ("T]) derivatives which lead to an axial 
velocity type multiplier which is always non- zero and relatively large for 
the cases to which this computational procedure is to be applied. 


Equation (71a) is then recast, using the appropriate transformations 
as discussed in the section on the computational plane, to obtain a differen- 
tial equation for h^ which is integrated at mesh points on the body during 
the computation. 


The metric coefficient h^ is then stored in an array and treated the same 
as any other flow variable. An exterior calculation then uses this data as 
input to trace streamlines on the surface of the vehicle, create the pseudo- 
stream surfaces (by taking noimals to the body at points along the streamlines 
of length 3 x flat plate boundary layer thickness), and interpolate in three 
dimensions to obtain all quantities and derivatives of interest in these 
surfaces. Normal derivatives are formed numerically. 


The method used for tracing streamlines on the body does not rely on simply 
stepping in the velocity direction (thus in most cases off the body) and then 
arbitrarily "pulling" the new point to a "corresponding" position on the body. 
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Instead, a method was devised to work directly with the angular position (0') 

s 

of the streamline on the body. A differential equation was set up for d0'/<iz 

s 

on the body, and was integrated in. a marching fashion. 

/s 

A /*. . 

Consider the vector triad (N,B^T) as shown in Figure 2o, N is the normal 
supplied analytically by our geometry package), B is a body tan- 

A 

gent lying entirely in the cross sectional plane, and T is also tangent to the 

A A 

body, and equals B x N. 

A ^ A 

N=Ni+Ni + Nk ■ 

X y^" 2 

A A A 

A 

B=Bi+Bd + Bk 
X y 2 

AAA 

T=Ti+Tj + Tk 
X y 2, 

A 

Since B must lie in the x, y plane, B h 0. Also, 

2i 

N*B = 0= BN + BN 
XX y y 


and since B is a unit vector. 


2 2 

1 = B + B 

X y 


Accordingly, 


= -V 


X y 


B = N / 

y x' 1 X y 


Finally, since T = B x N, 


T = B N 
X y z 


T = B N 

y X 2 


T = B N - N B 
2 X y X y 
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Now define the velocity directions (u, v, 'w) in the (N, B, T) directions. 

A 

Since v = ui + vj + wk 


= V*N = 0 = uN + vN + wN 
X y z 


V = V*B = uB + vB 

X y 


w=V.T=uT +vT +WT 
X y z 


In the (r,* 0,' p,') and (xj y,' z') coordinate systems (Pig. 5, h' = z), 
we define an elemental length (ds) lying in the x/ y’ (or x, y) plane, such 
that, 


We also note that. 


Accordingly, 


We can also write 




- dz' ", dz " 

W = -z-r- k ’ = -rr- k 

dt dt 


V _ ^ ^ ds 
w ~ dt dz ~ dz 

ds = v/w dz 


ds = dr'^ + r'2 d9'^ 


which yields 


and thus. 


de' = ds/W (dr’/de')^ + 


(v/w) 

(dr'/de)2 + r'^ 


Then, using a second order backwards integration, assuming (d0s/dz) 
linear between computational planes (see Figure 2j), 


''s ■ 

~ \ ^ )o [ (~^)l ' \~^ ) o 
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we obtain: 



This method of integration is used to trace any set of streamlines from 
specified positions on the body in the initial data plane. The radial posi- 
tion of these streamlines is, of course, supplied immediately by the geometry 

package for a given 0 ' . 

s 

Sonic Boom Data . Sonic boom calculations make use of the same data as the 
boundary layer quantity calculation. Another calculation performs a two- 
dimensional second order interpolation in given data planes to obtain the 
values of flow field variables (p, S, u, v, w) on a data cylinder enclosing 
the body, of user- specified radius, whose center line is the z or body axis 
(FRL) . 
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CONTROL POINTS 


Figure 24. - Exan^jle of modular breakup for aerodynamic coefficients. 



Figure 25. - Definition of coordinate system for metric coefficient 
calculation. 
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Figure 26. - Definition of coordinate system for streamline tracing technique 



Figure 27- - Diagram for streamline marching procedure. 








TYPICAL RESULTS 


During the course of this work the flow fields about a large number of 
geometries have been studied. Firstly to be sure that the numerical techniques 
used worked for all these configurations, second to compare the results 
with existing data, and finally to study flow field phenomenon which can't be 
studied through other means. In this section some of the results of these 
calculations will be discussed. Most of these results have been presented 
previously in references l8 and I 9 . 

In figure 28 the surface streamlines and pressure distribution are shown 
for an 80° slab delta wing at 30° angle of attack with = 9-6. At this 
high angle of attack the cross flow velocity expands around the wing tip and 
becomes supersonic on the lee side. This velocity component must vanish at 
the leeward symmetry plane and therefore a cross flow shock (Fig. 19a) is gen- 
erated. The strength of this shock at the body is demonstrated by the stream- 
line deflection across it. Figure 26 shows an axial and circumferential sur- 
face pressure distribution compared with experimental data (ref. I 8 ). The 
jump in pressure (circumferential distribution) at about 9 = 170° is due to 
the cross flow shock. The calculation of this flow from a starting plane 
(Z/Rj^ ^ . 8 ) to the end (Z/^ = i5) took about 15 minutes on the IBM 370/l68 
computer. 

Figures 29 , 30 , 31 and 32 all describe the flow about an early version 
of a shuttle orbiter configuration flying at an angle of 30° for M =26.1. 

CO 

The calculation was performed for an ideal gas, Y = 1.12. Figure 29 shows 
the top and side views of the shocks and body, while figure 30 shows the 
cross sectional views. The shock pattern in this calculation is quite com- 
plex, a strake shock is generated at Z ^ 38 O and intersects the bow shock 
at Z 440, a crossflow shock is generated starting at Z ^ 550 and a wing 
shock starts at Z ^ 8 OO and intersects the bow shock at Z ^ 85 O. Figure 31 
shows the surface pressure variation with surface distance around the vehicle 
at several cross sectional stations. Between Z = 8 OO and Z = 85 O there is 
a large pressure raise (Fig. 31 ) <iue to the beginning of the wing. The large 
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peak in pressure at Z = 85O could not be computed without concentrating grid 
points in the wing tip region. Near Z = 85O the wing shock intersects the 
bow shock and causes an expansion fan which reduces the pressure (at the wing 
tip) further downstream. This phenomenon can be observed up to Z = IO5O; 
while the peak pressure remains unchanged, the expansion moves toward the 
windward symmetry plane. The drop in the peak pressure between Z = 1050 and 
Z = 1100 is due to the wing tip turning parallel to the flow. Figure 32 
shows the radial entropy distribution at three circumferential positions 
(X is the noimalized radial coordinate). In these figures the thinning of 
the entropy layer is shown. This calculation (from Z = 50 ' to Z = 1280) 
took approximately 1.5 hours on the IBM 370/l68 using a maximum of 20 x 30 
grid points. 

Figure 33 shows the surface pressure distribution in the windward sym- 
metry plane on the forward portion of the O89-B shuttle orbiter. The calcu- 
latL ons were performed using equilibrium air thermodynamics (at 215 thousand 
feet altitude) and ideal gases at Y = 1.12 and 1.4. In each case the Mach 
n-umber was 26.1 and the angle of attack was 30° • The trends in the y = 1.12 
and real gas cases look very similar, while in the y = 1*^ calculation the re- 
compression (after the nose expansion) seems weaker. In the real gas calcula- 
tion the computer running time was increased by approximately 30 percent. Figure 
34 shows the surface streamlines on the O89-B shuttle orbiter. In the figure, 
the strength of the cross flow shock is demonstrated by its deflection of the 
surface streamlines. 

The windward symmetry plane surface, pressure distributions on a modified 
version of the current l40-C space shuttle orbiter are shown in figure 35 for 
M = 10.29, Y = 1-^j s^nd a = 20 and 25°. The vehicle was modified for the 
computations by increasing the wing sweep from 45° to 55° (fig. 35) to avoid 
subsonic axial Mach numbers near the wing tip which would occur for this 
value of Y* This modification has little effect on the windward symmetry 
plane pressure distributions. The calculated results are compared with ex- 
perimental data (unpublished) obtained at NASA/ Ames. The agreement is very 
good except near the trailing edge. The under prediction near the trailing 
edge is due to a mismatch in lower surface slope (of approximately 5°) between 
the experimental and numerical geometry models. 
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The circumferential pressure distribution for this configuration at Z/L = 
0.3 is compared with experimental data in figure 36 . Again the agreement' is 
very good. These calculations required approximately one hour of computing 
time on the GDC 66 OO using a maximum of 15 x 32 grid points. 

Another complex shock pattern is shown in figure 37. The vehicle is a 
fighter aircraft and the flight conditions are M =2.5 and o' = 6 °. In this 

CO 

flow field there was a canopy shock, a wing shock and an additional shock due 
to a recompression after the canopy. Figure 37c shows a series of cross sec- 
tional views in which the intersection of the canopy and bow shocks is evident. 
This calculation demonstrated that flow fields containing multiple imbedded 
shocks in a cross section can be computed treating all shocks explicitly. 

Figure 38 compares the computed and experimental surface pressures on 
this vehicle (M^ = 2.2, ' or = 5 and 10°). The experiment was run by Gmjmman 
Aerospace Corp. at KASA/Ames and the data are unpublished. This calculation 
(from Z = 0 to Z = J 45 ) took about 1 hour on the IBM 37 O /165 computer with a 
2h X 29 grid in each cross section. 

Figure 39 shows the inlet flow field for a supersonic fighter configura- 
tion. The axial station shown corresponds to the inlet forward lip station, 
so that this is the flow field ingested by the inlet. The maximum differ- 
ence between calculation and experiment is less than 3 percent. The experi- 
mental data were obtained from reference 31 and the calculation took approxi- 
mately 30 minutes on the IBM 370/165 with a 25 x 30 grid. 

Figure 40 shows surface pressure distributions on the X-I 5 aircraft top 
and bottom symmetry planes. In figure 40a the flight Mach number is 6 and 
in figure 40bthe Mach number is 4. The comparison with the experimental data 
is good. The deviation from the experimental results near the vehicle nose 
is due to the starting solution which was used (i.e., a conical flow solution 
was used to start the calculation). 

Figure 4l shows a sample of the type of sharp leading edge wing configura- 
tion which can be computed. The figure shows the computed shock pattern with 
the bow shock intersecting the wing shock. 
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Shown in figure 42 is another complex shock pattern. The vehicle is a 
hypersonic research aircraft (HSRA) configuration and the flight conditions 
are = 6 , Y = 1.2, and a = 0°. As shown in the figure the canopy shock is 
generated ahead of Z = 36 and is intersecting the bow shock at Z = 60 . A 
shock wave is generated by the vertical tail ahead of Z = 50 . Figure 43 
shows a lift vs angle of attack curve for the HSRA vehicle. Both the ex- 
perimental and Newtonian flow results were supplied by Mr. Lewis Clark of 
NASA/LRC. The figure shows that both the Newtonian calculation (with viscous 
effects included) and the present one compute lift accurately. The calcula- 
tion of the HSRA flow field (M = 6 a = 0°) took 1 hour on the IBM 3?0/l68 

CO 

computer using a 25 x 30 grid. 

The metric coefficients of the coordinate system based on the body 
streamlines (see figure 44) are one of the less, straightforward by-products of 
this computation. 

On normals to the body, flow quantities and their normal derivatives are 
also automatically computed (Fig. 45) and turn out to be smooth and satis- 
factory. 


/ 
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80° SLAB DELTA WING ATIVIoo=9.6, a = 30°&7 = 1.4 


STREAMLINES 



(a) 



SIDE VIEW 



Figure 29. - O89-B shuttle orhiter shock patterns. 





Z = 50.45 Z= 100.30 


Z = 150.16 



Z = 200.04 




Figure 30. - O89-B orbiter cross sectional shock patterns 
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Z = 550.74 


-72- 





Figure 31- - O 89 -B orb iter surface pressure. 
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Figure 36. - Circumferential pressure distribution, 
modified l4o-c shuttle orb iter. 
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M«,= 6„ 


7= 1.4 


Moo = 4., 


7= 1.4 



X-15 e = LENGTH OF VEHICLE (b) 


Figure UO. - X-15 surface pressure. 



Figure 4l. - Sharp leading edge wing configuration shock pattern. 
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CONCLUSIONS 


1. A computer tool has been developed for a high-quality computation of 
the inviscid flow field about realistic vehicle configurations at 
supersonic and hypersonic speeds. 

2. The underlying numerical techniques, while available for simple bodies, 
have been expanded, adapted and integrated into a tool. For high- 
quality results, all shocks are treated as shock discontinuities and 
entropy layers are handled in a special way. Moreover, a crucial element 
of the numerical technique evolved is the use of mappings to simplify 
the cross sectional shapes in the computational space and efficiently 
distribute a relatively small number of mesh points. 

3 . Equilibrium and frozen air are included. 

4. The conputer tool automatically produces specialized outputs such as 
aerodynamic coefficients, metric factors, normal derivatives in pseudo 
stream surfaces and sonic boom data. 

5 . The initial results obtained from the computer tool for orbiter and 
other configurations are quite satisfactory. 

6. As an orientation, the machine time expenditure for the flow field past 
a t 3 pical orbiter configuration at M = 26.1, a = 30°, and for a perfect 
gas is about I .5 hours on the C.D.C. 66 OO. 
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APPENDIX A 


Second Derivatives of Mappings 
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ITow: 


G - = W (W- )(W. ) + (W )(W )(W w. ) 

C/ 3vfg 3>72p Gg ? 4 

+ (W. )(W )(w ) 

5w4p IC 3 w4 

= - («c? - “cc 
“xp “ "" 

^ (Cop) 

Then we can get r r 6 „ and 6 ^ by taking the derivatives of 
xp’ yp’ x2 yp 

equations (7). In order to get r and 0 „ we take the 2 derivatives of equations 

zy. 

(8) and calculate x-_ and y_„ (again = A , = B , etc. since z = 2), 

2P PP PP zz PP , zz’ 

and we get from ( 9 ) 

* Vpp * --xp =‘p - ^x=^pp) 

■ - ('’y ^p * Vpp ^ ^xp^-p - ®x>^?p> 
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APPENDIX B 


Coefficients of Transformed Euler Equations 
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